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This paper considers the Monte Carlo dynamics of random 
dimer coverings of the square lattice, which can be mapped 
to a rough interface model. Two kinds of slow modes are 
identified, associated respectively with long-wavelength fluc- 
tuations of the interface height, and with slow drift (in time) 
of the system- wide mean height. Within a continuum theory, 
the longest relaxation time for either kind of mode scales as 
the system size N. For the real, discrete model, an exact 
lower bound of 0{N) is placed on the relaxation time, using 
variational eigenfunctions corresponding to the two kinds of 
continuum modes. 



I. INTRODUCTION 

Statistical models in which each microstate maps to an 
interface - which I will call "height models" - tend to be 
interesting, for the interface is often rough, and it turns 
out the model has critical correlations, jfj-^- This paper 
considers the corresponding dynamics., using one of the 
simplest lattice statistical models: the dimer covering of 
the square lattice [^,0. The dimers are placed on the 
bonds and every site must be touched by exactly one 
dimer. (This is equivalent, of course, to the packings of 
"dominoes" of size 2x1 furthermore the ground 

states of the fully frustrated Ising model on the square 
lattice llj] map 2-to-l to the dimer coverings. In this 
paper, I only consider the case where each packing has 
equal weight. The system is taken to be a square oi Lx L 
sites with periodic boundary conditions for the dimers. 

The statics can be described by mappings to free 
fermions, and exact solutions are possoble using Pfafh- 
ans 1^,0 or transfer matrices (l^Jl^ . There is a nonzero 
ground state entropy of 0.2916 per site d). The cor- 
relation functions are critical | [T5[ | (power-law decaying). 
This is easiest understood after mapping the dimer pack- 
ings (configuration-by-configuration) to configurations of 
"heights" z{r) living on dual lattice sites, r = {x,y), rep- 
resenting a rough interface in an 3-dimensional abstract 
space 0-^]. 

When such a model is endowed with dynamics, a cen- 
tral question is how the equilibration time t(L) (to be 
defined shortly) scales as a function of the system diam- 
eter L? In particular, what is the dynamic exponent z 
p6[ in t{L) ^ L'^l Since the static continuum model has 
a gradient-squared free energy, one would suspect that 
z = 2. Indeed, scaling with z = 2 has been seen numeri- 
cally in simulations of the square lattice dimer model p8| . 



(It has also been seen in antiferromagnetic Ising models 
and random-tiling quasicrystal models with height 
representations.) However, it has only been proven pc| ] 
that (with our definition of time scale) t(L) < 0{N^), 
which implies z <Q. 

This paper has two aims: (i) an explicit approximate 
description of the slowest eigenmodes of the time evolu- 
tion, based on continuum theory - this theory is appli- 
cable, with small changes, to any height model; (ii) the 
outline of a proof that t{L) > 0{L^) (and hence z > 2.) 
The 

The outline of the paper is as follows. The remainder 
of this section specifies the model, in particular the height 
mapping and dimer-fiip stochastic dynamics, and also 
points out an exact correspondence of the entire spec- 
trum of the quantum dimer model to that of this dynam- 
ics. I then take up a continuum model in terms of the 
height representation (Sees. ||) which is defined by stan- 
dard Langevin dynamics, verifying that this is the appro- 
priate coarse-graining of the microscopies. This produces 
an approximate description of all the slow eigenmodes 
(Sees. H and HI). Intriguingly, the slowest mode in any 
finite-size system is not the longest-wavelength capillary 
wave, but a "height-shift" mode involving diffusion of the 
system- wide average "height" direction (Sec. HI). 

Finally (in Sec. ^ I explain how rigorous bounds on 
the dynamics can be proven. The fundamental concepts 
for this proof are (i) the use of Fourier analysis to par- 
tially diagonalize the evolution matrix W, (ii) construct- 
ing variational "wavefunctions" guided by the results of 
the earlier sections, and (iii) taking advantage of the 
Fourier spectrum of height fluctuations derivable from 
the exact solution of the model. In the conclusion, I 
mention other spin models to which these results should 
be relevant. 



A. Height representation 

The explicit rule for constructing a height pattern 
{z(r)}, given a snapshot of the dimers, is shown in Fig. 
1: z(r) — z(r') = —3 if there is a dimer between r and 
r', or -1-1 if there is no dimer, where the step from r to r' 
is taken in a counterclockwise sense about the even sites. 
The net height difference is zero for the path around one 
plaquette, and hence by induction for any closed path, 
showing the consistency of the definition. However, it 
is nonunique in that adding the same constant to each 
^(r) makes an equally valid height representation of the 
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same dimer configuration; I shall fix this constant after 
defining the dynamics (below). 

It is necessary in writing the continuum model for 
z{r) (e.g. to define its Fourier transform) that z(r) sat- 
isfy periodic boundary conditions; but periodic boundary 
conditions for dimers only imply z(L,y) — z(0,y) = Wx 
and z(x,L) — z(x,0) — Wy, where the "winding num- 
bers" Wx and Wy are multiples of 4, so that the sys- 
tem has a mean tilt is {wx/L^Wy/L). Local update 
rules (such as I am about to define) conserve the wind- 
ing numbers. Thus one could define substracted heights 
z'(x, y) = z(x, y) — [wx/ L)x — {wy/L)y which do satisfy 
periodic boundary conditions. In the limit of small tilts, 
z'(r) would obey the same continuum dynamics.]^ This 
trivial generalization is not worth the added complica- 
tion in notation; in the remaining sections I will only 
configurations with {wx,Wy) — (0,0) 

The ground states of the fully-frustrated Ising model 
on the square lattice may be mapped (2-to-l) to dimer 
coverings of the dual square lattice. The rule is: sim- 
ply draw a dimer across every violated bond. Therefore 
that model has a height representation and all my results 
apply equally to the fully- frustrated Ising model. 

The height variable will play the role of the "order 
parameter" in this paper. The more customary order 
parameter, for any of the height models, is a spin opera- 
tor m{r) (or a dimer occupation operator, in the present 
case). However, as outlined in Appendix ^), such an 
m{r) is simply a sinusoidal function of the local height. 
Since the long-distance and long-time height fluctuations 
discussed here are Gaussian, it is possible from them to 
compute the correlations of m(r) which turn out to be 
algebraic with nonuniversal exponents. 

B. Dynamics 

It is possible to turn any dimer covering (of zero mean 
tilt) into any other one by a succession of "dimer flips" 
each affecting two dimers on opposite edges of one pla- 
quette (see Fig. 1(b)). Thus the model is endowed with 
a stochastic (Monte Carlo) dynamics in continuous time 
as follows: select plaquettes at random, at a rate N per 
unit time where N = L'^ is the number of sites (i.e. on 
average visit each plaquette once per unit time). Flip the 
plaquette if it has two dimers as in Fig. 1(b), otherwise 
do nothing. 



On the other hand, the natural zero-temperature dy- 
namics of the fully frustrated Ising model is to choose 
a spin at random and flip it if the energy change would 
be zero, otherwise do nothing. This induces exactly the 
same dynamics on the dimer configurations as specified 
in the preceding paragraph. 

To eliminate the arbitrariness in defining z(r) and en- 
sure that the coarse-grained dynamics is continuous in 
time, we relate z(r,t) at different times, by specifying 
that a dimer fiip on a plaquette changes only the z{r) 
value in that plaquette's center. (Notice that the possi- 
ble z(r) values on a particular site can only differ from 
the initial value by a multiple of 4; the values of z{r) 
(mod 4) define four fixed, square sublattices. |^). 

Now let {pa{t)} be the probability of being in mi- 
crostate a at time t. The master equation states 

^ = E iPf^ - p^) -T.^^m (1) 

where (a —f [3) means summing over configurations (3 
which differ from a by one dimer fiip; thus 

^ FaSct3 - Aa[3 (2) 

where A is the adjacency matrix (elements unity if a 
and P are related by a fiip, and zero otherwise), and 
Fa = J2f3 Aaf3 is the number of "fiippable" plaquettes in 
configuration a. These matrices are Mq x Afo, where Mq 
is the number of microstates. It is well known that the 
matrix W has nonnegative eigenvalues (via the Perron- 
Frobenius theorem, since W is a stochastic matrix). The 
eigenvector of zero eigenvalue is Pa = 1/Mq, the weight 
of the (equilibrium) steady state (which is unique, since 
dimer fiips connect all microstates). Furthermore, any 
time correlation function in the system can be resolved 
into a sum over eigenvalues A of W, in the form 

EcAe-^' (3) 

A 

Then it is normal to define the system's equilibration 
time 

r{L) = Kmn~\ (4) 

where X'^^^ is the smallest nonzero relaxation rate. The 
equilibration time t{L) of the system is defined as the 
inverse of the second smallest eigenvalue of W. 



^ If tilts are nonvanishing as L — * oo, the free-energy func- 
tional (^ must be generalized to have different stiffnesses for 
components of the gradient parallel to and transverse to the 
tilt direction; the same power laws for the dynamics would be 
deduced. 



C. Quantum dimer model 

In the quantum dimer model, the basis states are taken 
to be the dimer coverings, and the Hamiltonian is taken 
to have matrix elements 

Hap = -Acf} + VFaSap. (5) 
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The first term describes dimer flips ~ like a particle quan- 
tum "hopping" with amplitude t on the microstate graph; 
the second term is a "potential energy" which penalizes 
each flippable plaquette. 
When V — 1, obviously 



K 



a(3 



(6) 



Ref. noted that implies the ground state wave- 
function of (|^) is a superposition of all the dimer packings 
with equal amplitudes. In other words, that (quantum) 
wavefunction is proportional to the steady-state proba- 
bilities of the (classical) dynamics (Q) . 

Here I note that (^) further implies a one-to-one corre- 
spondence of all the eigenstates of (|^) to normal modes of 
the master equation. Thus the bounds derived in this pa- 
per for the slowest relaxation rate are equally valid for the 
energy gap in the quantum model, and the approximate 
eigenf unctions and eigenvalues found in Sees. || and [II 
also describe the low-energy spectrum of the quantum 
dimer model. 



II. CONTINUUM THEORY OF DYNAMICS 

This section develops the continuum ("coarse- 
grained" ) version of the dynamics. In this form the model 
has an easily visualized physical meaning and (being lin- 
ear) is solvable by standard and almost trivial techniques. 
In general, the slowest modes are associated with relax- 
ation of the "order parameter" 1 16 2^ ; in the dimer cov- 
ering model, the height variable plays the role of a (hid- 
den) order parameter, so the dynamics are phrased in 
terms of it. 

It should be noted that this theory is general to all 
height models^ The only specific information from the 
dimer model is the numerical value of the elastic constant 
K and of the height space lattice constant , which are 
also known for many other height models. 



A. Continuum equations and Fourier modes 

First I review the coarse-grained picture of the height 
dynamics. The static free energy functional has the form 



F 



d\—\Vh{r)\^ 



(7) 



Here h{r) represents a smoothed version of z(r) and K 
is the stiffness constant controlling the fluctuations of 
the "interface". (From here on I assume zero net tilt 



^ In the cases that 2:(r) has more than one component, it is 
necessary to let F be a tensor. 



of z{r) and of h{r) so these variables satisfy periodic 
boundary conditions.) The fact that the dimer model is 
rough (described by ^ is nontrivial: several other height 
models, defined in similar ways, are found to be smooth 
in which the interface turns out to be smooth, 
or marginal. The roughness is confirmed only through 
the calculation in App. 

The customary dynamics for such a field theory (see, 
e.g., [p3|) is formulated as a Langevin equation. 



dhjr) 
dt 



Sh{r) 



C(r,t) 



(8) 



Here F is the kinetic (damping) constant measuring the 
linear-response to the force 6F/6h{r), and ^(r, t) is a ran- 
dom source of Gaussian noise, uncorrelated in space or 
in time. In order that Eq. (^) have for its steady state 
exp(— F) with F given by Eq. (|^), the usual condition 



(C(r,t)C(r',0) =2F<5(r-r')<5(t-0 



(9) 



must be satisfied. 

In the context of models of real (rough) interfaces of 
cry stals, Eq. (||) is known as the Edwards- Wilkinson 
process . There is a large literature on more elaborate 
equations of this form (usually with nonlinear terms) [ p4[ . 

To make plausible the assumption of uncorrelated 
noise, one must consider the action of the microscopic 
dynamics on the microscopic heights z{r). An elemen- 
tary dimer flip changes z{r) on just one plaquette by 
±Az = ±4, and the next dimer flip occurs at another 
random place. Thus the net height is not conserved, the 
change is local, and uncorrelated in time, which are mod- 
eled by the identical properties of the Langevin noise in 
(I). 

We can Fourier transform the above equations, since 
periodic boundary conditions maintain translational 
symmetry. My short-distance cutoff prescription for 
this continuum field theory shall be that the fields' 
Fourier transforms have support only in the Brillouin 
zone (— TT, 7r)^; in other words, the only allowed q values 
are those which are defined for the microscopic lattice 
model (see Sec. Then ^ becomes 



dhq 

~dr ^ 

with Gaussian noise 



(Cq(0*C;(i')>-2F<5q,q,5(t-t'). 



(10) 



(11) 



(Note that (5q,q' is the discrete (5-function, appropriate to 
the discrete lattice of wavevectors corresponding to pe- 
riodic boundary conditions.) Thus the different Fourier 
components are decoupled in (|l^ . For each of them (ex- 
cept for q = 0), ( p^ ) is a one-dimensional Langevin equa- 
tion with a restoring force. 
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B. Dynamic scaling 



probability density ^'({/iq(i)}):0 



The correlation function of the heights is easy to derive 
from the Langevin equation: it is 



1 



.-Ah(q)t 



where the relaxation rate is 



A„(q) = ri^|q|- 



(12) 



(13) 



Hohenberg and Halperin argued [Q, that the dy- 
namic exponent is best defined by the relation between 
relaxation rates and wavevectors. Then ( |l^ ) implies 
z = 2. As a corollary to (p^), the smallest relaxation rate 
of a Fourier mode corresponds to the smallest nonzero 
wavevector q^min = (27r/L,0) i.e. 



A,,(q™„) = (47r2ri^)i- 



(14) 



a plot of A/i(q) versus q revealed the behavior (|13|). A 
similar method was used in a study of the Ising an- 
tiferromagnet of general spin at T = on the triangular 
lattice, which also has a height representation. 

Ref. |l^ simulated the random dimer model with 
dimer-flip dynamics; however rather than periodic 
boundary conditions, that work used "Aztec diamond" 
boundary conditions |0 such that h{r) is fixed (and 
spatially nonconstant) along the edges. The continuum 
equations below would still predict r(i) ~ in that ge- 
ometry, as was observed in the simulations. ||l^ Note, 
though, that they define t{L) without the benefit of 
Fourier analysis, as the mean time it takes two (initially 
independent) replicas of the system to coincide, when 
evolved using identical random numbers sequences. 



C. Fokker-Planck mode spectrum 

Just as the discrete stochastic dynamics implies (|^), 
the continuum stochastic dynamics (^) implies the fa- 
miliar Fokker-Planck equation for the evolution of the 



d \ ^ 



dh, 



dh. 



K\cifh^)P{{h^{t)}) (15) 



Random tilings, with vertices not constrained to lie on 
a periodic lattice, are studied as models of quasicrystals 
pTf . These too can be mapped to effective interfaces 
z{y). (A complication, unimportant for the present dis- 
cussion, is that the height function z(r) or /i(r) has two 
or more components in the quasicrystal cases.) In one 
case of a quasicrystal random tiling (in three spatial di- 
mensions), (|l^ ) was confirmed by simulation []l9|] . After 
z(r, t) was constructed, the data were numerically Fourier 
transformed to give Zq{t) at selected (small) wavevectors 
q; for q small, that is essentially hq{t). Then the time 
correlations {zq{0)zq{t)) were fitted to the form ([13) and where 



A physicist analyzing the Edwards- Wilkinson dynam- 
ics (H) would usually have stopped at (|lj) . That does in- 
deed represent the slowest Fourier mode at each wavevec- 
tor, but there are many more eigenmodes of the general 
equation (|l5|). These modes are worth computing be- 
cause (i) they permit computation of more general cor- 
relation functions than ( p^ ) (ii) the analytic form of the 
modes, derived below, might inspire improvements on the 
variational eigenfunctions used in Sec. ^ and (iii) these 
modes correspond to excited states in the quantum dimer 
model at its critical point (see Subsec. ttO, above). 

The unique zero eigenvalue of Eq. (|15|), corresponds 
of course to to the Boltzmann distribution which is a 
Gaussian: 



(16) 



q#0 



To construct all the other eigenfunctions, it is conve- 
nient to write ^i{hq{t)}) = P{{hq{t)})Poi{hq{t)})-^/^ 
so that the time evolution operator of ^'({ft.q(i)}) is Her- 
mitian: 



dt 



(17) 



= ^r/^|q|2CtqCq 

where the "annihilation" operator is 



(18) 



Cq ^ (lif |qn-i/2( + iXlqp/iq) (19) 

and the corresponding "creation" operator is 



C\ ^ (ii^|qn-V^( _ + iK\q\'h^q) (20) 
dhq 

Of course C^q commutes with Cq' for all q' ^ q. Ob- 
viously this is mathematically identical to the quantum 



^ The customary abuse of notation is committed in which 
hq appears to be manipulated as if it were a real variable; the 
convention dh* /dh — justifies the manipulations in ( p^ ) - 

(§)• 
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Hamiltonian for a set of harmonic oscillators with fre- 
quencies A;i(q) given by ([l^), with a "ground state wave- 
function" = Pf 



1/2 



Now we can write any other eigenfunction: 



(21) 



where {n(q)} are any nonnegative integers, correspond- 
ing to the occupation numbers of the oscillators. When 
translated back in terms of P{{hq{t)}), we see that each 
eigenstate is a product of PQ({hq(t)}) , times polynomials 
in {(iif |qp)^/^/iq}. For the "elementary excitation", in 
which n(q) = 1 for one wavevector and zero for all the 
others), this polynomial is exactly hq; this explains why 
the correlation function ( |l2| ) sees only that one eigen- 
mode. 

The eigenvalue corresponding to (|2l|) is 



(22) 



Ato^({?^(q)}) = ^7i(q)Aft(q) 



corresponding to the total energy of the quantum oscil- 
lators. The net wavevector is 



E 

q#0 



n(q)q. 



(23) 



Notice the many degeneracies resulting from the fact that 
A(q) = A(— q): not merely the degeneracies due to global 
symmetries such as q — > — q, but less trivial degeneracies 
such as the one between the mode with n(q) — n(— q) = 
1 and the one with n(q) — 2, n(— q) = 0. 



Gaussian random walk. Q 

When h{r) describes a genuine interface, states with 
different h are all distinct, and the distribution of h sim- 
ply spreads diffusively without ever reaching a steady 
state. However, in the dimer model (and all other "height 
models" [0j5] ||]) the height map is one-to-many: here a 
global shift of z{r) z(r) + 4 describes exactly the same 
dimer configuration, so the image space of h{r) should 
be considered a circle of diameter 4. Thus the distribu- 
tion function, which begins sharply peaked at a particular 
value of h, will evolve to a uniform distribution at some 
rate. 

To make this mode more concrete, it may help to com- 
pare with the behavior of a spin operator as seen in 
App. A system with a height distribution peaked at 
(say) 1 and 3 has more dimers in one orientation than 
in the other. As this n(0) = 2 mode decays, this polar- 
ization of orientations will decay; thus these modes have 
quite real physical meanings. 

The random walk behavior (analog of (l2h) is 



{\hit)-hm')^D{N)t 



with a diffusion constant 
D{N) 



2r 



(25) 



(26) 



It is interesting that although the continuum theory can- 
not provide the numerical value of the coefficient F, it 
does predict an exact ratio (in the limit N — L"^ ^ oo) 
between the relaxation rate of the modulation h{qmim t), 
and the diffusion rate of the wandering of h{t). 



III. HEIGHT-SHIFT MODE 

Besides the translations in real space, a height model 
(as I define it) has the additional symmetry of trans- 
lations in "height space" (the target space of z(r) and 
h{r)). Correspondingly we will find another kind of slow 
mode (in a finite system), corresponding to a random 
walk of the mean height, which I will call the "height- 
shift mode" . It will turn out to be the slowest mode. 

Consider the average height 



h{t) EE 



(fr h{r,t) 



(24) 



i.e. N'~^^'^ho{t) (the factor N^^^ comes from my normal- 
ization convention for Fourier transforms.) Of course, 
this is just the q = mode that was excluded from in the 
preceding section (e.g. in ( |2l| ) and (p^)). The Langevin 
equation (0) with |q| = tells us /io(0 simply executes a 



A. Fokker-Planck modes of h{t) 

The Fokker-Planck equation for h is simply the diffu- 
sion equation; its eigenfunctions are simply plane waves 
as a function of h, 



ipqih) = ex.p{iQh) 
and the corresponding eigenvalue is 



(27) 



(28) 



with D{N) given by (p6|). If h were diffusing on a line, 
then any Q would be valid in (28) giving a continuum of 
eigenvalues. 



* This behavior (and the assertions deduced from it) are 
valid only while the height model is in a "rough" phase. One 
motivation for understanding the relaxation modes (or quan- 
tum energies) is that they might serve as a diagnostic to dis- 
tinguish rough and smooth phases in simulation results. 
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But, as noted above, a global shift z{r) z{r) + 4 of 
a microstate corresponds to the same microstate. Thus, 
the only modes which can correspond to modes in the mi- 
croscopic model are those periodic under h{r) /i(r)-|-4; 



g-n(O)- (29) 
for any integer n(0). and the smallest such eigenvalue is 

(30) 

Thus the complete set of eigenfunctions are 



TT^r 

A7-(vr/2) ^ 



^'({/iq}; {n(q), q ^ 0})V'„(o)V2(/i) (31) 

a product of (|2^) and (p7|), with any integer n(0). The 
corresponding eigenvalue is 



Atot(Wq)}, q ^ 0) + n(0)2A^(7r/2) 



(32) 



with the first term from (p^). The eigenvalue ( |30| ) given 
by n(0) = ±1 and 7T,(q) = otherwise, is in fact the 
overall smallest nonzero eigenvalue.^ Thus (^0|) is smaller 
than ( p^ by a factor (IGiiT)"^. From Appendix |^, the 
exact value is if = 7r/16 so this ratio is I/tt. 

Is the height-shift mode qualitatively distinct from the 
capillary wave modes? Clearly there is a close connection 
- the uniform fluctuation of the entire sample of diam- 
eter L is much like the fluctuation of one quadrant of a 
system of diameter 2L with a long wavelength Fourier 
mode. I would argue there is a distinction, in that the 
rate of the hit) relaxation depend on "internal" details 
of the height model; it cannot be inferred given only the 
capillary spectrum. 

For example, I have noted the fully-frustrated Ising 
model with single-spin-flip dynamics is identical to the 
dimer model except that it has a height space period of 8 
rather than 4. Thus one must replace 7r/2 7r/4 in ( |29| ) 
and the smallest eigenvalue (^^ is smaller by a factor 
1/4. (In the fully-frustrated Ising model, a shift z z+4 
reverses all the spins; thus the mode with Q = ±tt/2 is 
the slowest mode that is odd in spin.) In the other direc- 
tion, I cannot rule out the possibility that in some height 
model, the stiffness constant K and height periodicity 
might have numerical values such that Kaf^ > 1, in which 
case the slowest mode would be the Fourier mode. 



The relationship of the height-shift mode to the cap- 
illary wave modes (both of which are a consequence of 
height-shift sysmmetry), is reminiscent of the relation- 
ship between two kinds of low energy excitation in a 
quantum spin system (both consequences of spin rota- 
tional symmetry.) The height-shift mode is the ana- 
log of the change in total spin number (giving energies 
5(5' -I- l)/2x); the capillary mode is the analog of a spin- 
wave mode. 



IV. CORRESPONDENCE TO DISCRETE MODEL 

The low-lying modes of the discrete model should 
be well described by the "quantum numbers" n{c{) of 
the capillary modes and height-shift modes of contin- 
uum model. Now, Po{{hq{t)}) corresponds to the 
eigenmode of the discrete model which has equal weight 
Pa "^ = I/Mq for every microstate. Thus the prescrip- 
tion for constructing the approximate eigenmode of the 
discrete system is 



J'a;n(q) 



^({^q};{''^(q)})V'«(o)V2(^ 

*0({2q}) 



(33) 



Here Zq is considered an implicit function of the discrete 
state label a. We obtained (^^ simply by replacing 
Zq, and /i — > z defined by 



(34) 



which one expects to be valid for small enough n(q) and 

q 

However, ( p3| ) does not map every mode of the con- 
tinuum equation to an approximate mode of the discrete 
one. Indeed, since each n{q) is unbounded from above, 
the continuum kernel Wh has infinitely many indepen- 
dent eigenmodes, whereas the discrete kernel W has only 
Mo eigenmodes. Presumably, when ri(q) is so large that 
n(q)|qp is of order unity for some q the anharmonic 
terms that I omitted in writing (^ become important 
and mix this mode with others. Thus only qtot is a good 
"quantum number" for labeling the higher modes. 

The low-energy eigenstates of the quantum dimer 
model (^) can be predicted from the above analysis. The 
wavefunctions should be approximated by (|3^) and the 
energies given by (^). These predictions could be com- 
pared [ p5| with the results of recent exact-diagonalization 
studies on 8 X 8 lattices p6[). 



^ Of course, ( po| ) would not be smallest in the case of h{r) 
fixed along the boundaries, as in and [^, since that 
forbids /i- wandering modes. Notice also that, for systems of 
unequal length sides, the ratio of (^|) to ( |30| ) is decreased by 
a factor Lmin/Lmax SO that for a sufficiently elongated system 
the Fourier mode is the slowest one. 



® The discrete model has additional lattice symmetries (rota- 
tions by 7r/2 and reflections), which correspond to additional 
"quantum numbers" ; it furthermore has a "locking" tendency 
(to favor a particular value of h{r) (modulo 1). Discussion of 
these complications will be deferred to a later publication [E5| . 
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V. MICROSCOPIC THEORY 

Now we return to the exact microscopic dynamics, as 
introduced in Sec. I. The microstates are viewed as nodes 
of a graph in an abstract space, with each possible dimer- 
flip (and its inverse) forming an edge of the graph: the dy- 
namics described by (|]) is a random walk on this graph. 
As noted already, the graph of microstates with tilt zero 
has only one connected component. The matrix W is not 
only stochastic but symmetric. 



A. Variational bound 

The key idea in this section is the use of a varia- 
tional guess for the eigenfunction. This is mathematically 
equivalent to the variational bound on t he e igenenergy in 
the quantum dimer model (see Subsec. 1 C.) 



For any vector {(j)a}, the smallest eigenvalue A„ 
isfies 



A. 



<Xd> = 



sat- 



(35) 



Here (/>* means the hermitian conjugate. From (H), the 
numerator is 



(36) 



{af3) 



where (a/?) means every pair of microstates, connected 
by a spin flip, is counted once. 

In the two cases considered in this paper (following 
subsections), 1 0a — 0/3 1 = |A0| turns out to be the 
same for every flip move. Then the sum (|3^) reduces 
to A0p Fa (the 1/2 cancels the double-counting of 
each graph edge), and finally to ^fN\A(j)\'^ Here / is the 
probability that a given plaquette is "flippable" , i.e. fN 
is the average "coordination number" of the microstate 
graph. The exact value is 



/ - 1/8. 

Meanwhile, the denominator in ( ^5|) is just Mq{ 
Thus finally the upper bound is 



(37) 



A 







1 |A0P 
2-' (\m 



(38) 



When applied to the entire set of eigenvalues of W, 
eq. (^8|) is not very interesting, since we already know 
that Xmin = (the eigenvalue of the steady state) . How- 
ever, the variational argument (and every equation in this 
subsection) is also valid when restricted to a subspace or- 
thogonal to that of the ground state. Then, the bound 
(p8|) can be useful since the Xmin of such a subspace is 



usually positive. Q Indeed, the overall smallest nonzero 
eigenvalue - the goal of this paper (see (||)) - is expected 
to be the Xmin of one of these subspaces. 

The above variational argument was first presented p^ ] 
with a physical interpretation in terms of the normalized 
dynamic correlation function 



C4.it) 



(39) 



Then (7^(0 ^ 1 - i(|0(O) - (j){t)\^) / {\(f>\^); but at short 
times, (10(0) - 0(t)P) = {\Acj)\'^)fNt, since there are on 
average fN independent places where a flip could occur. 
Thus the upper bound on Xmin can be rewritten as the 
initial decay rate of the correlation function, ||27[] 



A. 



dt 



(40) 



Ref. (and similarly |2g]) applied (^0|) to bounding 
the dynamic critical exponent by a function of the static 
exponents (for a system with nontrivial exponents) . 

In some circumstances, A<^ could be a good estimate 
of Xmin- It would bc exact if subsequent steps in were 
uncorrelated. (We assume our variational 0q. always has 
nonzero projection onto the slowest mode.) Then 0(t) 
performs a random walk, and C^^it) = exp(— A^i), a pure 
exponential decay. In light of (|^), in which coefficients 
allowed by symmetry are expected to be generically pos- 



itive, we would obtain A,, 



A. 



However, dynamics actually adopted (in eq. (|l|) is such 
that the steps are manifestly anticorrelated in time. For, 
if a plaquette flips, it must undergo the reverse flip the 
next time it flips, unless its configuration has been shuf- 
fled in the meantime by flips of the adjoining plaquettes. 
So one would expect (but did not prove!) that Xmin < A^ 
as a strict inequality. 



B. Fourier modes 

First we review the use of basis states with definite q 
vectors. Consider the action of a translation of r; this 
induces a permutation of microstates which manifestly 
commutes with W. Thus all eigenvectors must have a 
definite wavevector q (or can be chosen thus); that is, 
a translation by r simply induces a multiplication by 
exp(iq • r) of the eigenvector. Furthermore W must have 
zero matrix elements between vectors with different q; 
thus W becomes block-diagonal in the new basis of states 



^In subsequent subsections, labels will be attached to Xmin 
to indicate which subspace they belong with, but the labels 
have been omitted in the general argument here. 
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with definite q. Hence ( p35| ) is vahd for the smallest eigen- 
value Amm(q) in each block of definite q; that is, if {0q} 
has wavevector q, then 



(41) 



A suitable variational state-space vector (for q 7^ 0) is 
suggested by th e "elementary excitation" eigenfunctions 
of Subsec. II C ; they consisted of the steady-state eigen- 
function ^'o multiplied by hq. The microscopic analog 
of the gaussian is the equal-weighted eigenfunction of 
yV, hence we choose 



(42) 



(Here Zq{a) means take the configuration z(r) corre- 
sponding to microstate a and Fourier transform it for 
wavevector q). 

If a, are related by one dimer flip on the plaquette at 
rf, the corresponding configurations of z(r) differ only by 
±Az = ±4 at r/. Consequently |A0P = \Az\^/N. On 
the other hand, the denominator of (|3g) is exactly 

^|zq(a)p=Afo(|z(q)P) (43) 

a 

(recall each microstate is weighted equally). Thus the 
variational bound has the form 

A™„(q) < 8//(|z(q)p) (44) 

Finally, we know (|z(q)p) = 1/-R r|qp at small q, as de- 
rived in Appendix |b[ Eq. ( B10| ) . This gives the result 
(for small q) 



A™„(q) < 8fK\q\- 



(45) 



If we assume 
coefficient: 



, then (p5|) gives a bound on the kinetic 



8/ 



(46) 



This approach - constructing a variational vector from 
the Fourier transform of a local operator - is equivalent 
to the "single-mode approximation" used for quantum 
many-body systems such as superfluid helium p9| ]. In- 
deed, already noted (in the context of the quantum 
dimer model) that an exact upper bound on the eigen- 
value is implied. Their proposed variational wavefunc- 
tion is the Fourier transform of an dimer density operator 
nr{r). In fact, this is actually just ATz(r) (difference op- 
erator in the direction t = x,y)] thus their choice differs 
from (^ ) only by a derivative. They did not identify the 
excitations as capillary waves, but derived the exponent 
z = 2 using the correlations of pSl. n 



C. Height-shift mode 

The variational trick also works for the random walk 
behavior of (|2^). Obviously the microscopic analog of 
( p^ ) is 'z{t) defined by (jsj). Then the q = modes can 
be further block-diagonalized into modes with particular 



wavevectors Q in height space, as defined in Subsec. Ill A 



hence the variational bound is valid within each such 
block (to be labeled "(0;g)"). 

In analogy to (E^), we use the variational vector 



(47) 



In this case, |A0| = | sm{QAz/2N) \ for any state a, with 
Az — ±4 the same as before, and obviously {(pal — 
1 in any state. So from ( ^8| ) we obtain the bound 
A™„(0;Q) < i/7Vsin2(4Q/iV), or 



(0; Q) < 8/(7r/2) 



(48) 



for large N. 

Just like the results ([l3| ) and (|28| ) for the two respec- 
tive kinds of mode in the continuum model, the best 
variational upper bounds for slow relaxation rates in 
the microscopic model correspond respectively to |q| = 
|q„„„| = 27r/L in (||) or to Q = 7r/2 in (||). Either kind 
of bound implies a lower bound on the relaxation time r 
which is 0{L^). The better bound comes from (39), and 
after substituting the value of / (p^) gives 



T > (2/7r)2iV. 



(49) 



It is interesting to note that ([48|), when compared to 
( ^ ) and (28), gives the same bound on F as Eq. (|4^ ) 
derived from Fourier modes. This might be taken as an 
approximation for F, which amounts (as already noted) 
to neglecting the anticorrelation of fiips. Such an approx- 
imation for F was made in a random tiling quasicrystal 
model. (See Sec. IV A and footnote 36 of Ref. 0). 
This turned out ||l^ to overestimate the true F (esti- 
mated from the simulation) by about 50%. 



VI. DISCUSSION 
A. Summary 

The results in Sees. ^ and ^ exemplify the fruitfulness 
of Fourier analysis in models with translational symme- 
try; indeed two types of Fourier transform were used. 



* Caution: the correlation function exponent entering this 



calculation is always exactly 2, even when the spin operator 
exponents decay by exponents rj ^2 (see Appendix^. That 
situation arises in the other height models, or in the dimer 
model when confugurations are weighted unequally. 
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For capillary wave modes (Sec. ||, and subsec. |VB| ), 
the Fourier transform operates with respect to transla- 
tions in physical space. For height-shift modes (Sec. Ill 
and Subsec. VC), it operates on translations in "height 
space" . For each type of mode, two arguments have been 
given, one based on a coarse-grained field theory (Sees. i| 
and III) and the other based on exact bounds (Sec. 
which indicate the the longest relaxation time scales with 
system size as i^. 

The only rigorous consequence of the present argument 
is a lower bound (^9|) on the longest relaxation time. For 
computational purposes one prefers, of course, an up- 
per bound on the time needed to equilibrate the system. 
Towards this end, the present calculation merely sug- 
gests the possible usefulness of Fourier analysis in such a 
demonstration, and warns that the slowest mode is not 
a Fouri er mo dulation at all but the "z diffusion" (see 
Subsec 



III AD 



B. Computational physics 

This study is relevant to the dynamics of other spin 
models. In particular, its results apply directly to the 
critical dynamics of the fully-frustrated Ising model on 
the square lattice since that model (at its critical tem- 
perature, Tc = 0) maps to the dimer covering The 
single-spin-flip dynamics in the spin model maps to the 
dimer-flip dynamics used here. 

The derivations could be trivially adapted to the dimer 
covering of the honeycomb lattice and hence to the crit- 
ical (Tc = 0) dynamics of the triangular Ising antiferro- 
magnet (equivalent to the dimer covering [Q). The be- 
havior argued here can be expected in other spin models 
which have height representation and single-spin update 
rules, e.g. the T — three-state Potts antiferromagnet 
on the square lattice (equivalent to the 6-vertex model 
and hence the BCSOS model). In those cases, how- 
ever, no exact results or rigorous bounds for the static 
fluctuations are available to replace Appendix ^ 

Understanding the local dynamics addressed here is 
also a prerequisite to addressing the dynamics of height 
models endowed with loop-update or cluster-update 
rules, including fully-frustrated Ising models the 
three-state Potts antiferromagnet on the square lattice 
H, the BCSOS model @, and coloring models ||.[| 
More speculatively, it would be interesting to check 
whether there is any connection between the Swendsen- 
Wang |33| algorithm and the height representation in the 



case of the ferromagnetic 4-state Potts model, the parti- 
tion function of which can be mapped to that of a height 
model |36|. 

Incidentally, the static correlation exponents of height 
models have been obtained from Monte Carlo simula- 
tions much more accurately via Fourier analysis of the 
heights, than via direct measurement. It seems likely that 
the same is true for dynamic measurements; this might 
help in evaluating the dynamic exponent for cluster- 
acceleration algorithms (see Rcf. pO[) 

Somewhat analogous to the loop-update rules are the 
"zipper moves" [ p^ , ^ of certain random tiling quasicrys- 
tal models, such as the (two-dimensional) square-triangle 
random tiling. In that case, the correlation time was 
measured to be of 0(1) update move per site how- 
ever each update move involved 0{L^) sites, so the scal- 
ing of the net correlation time was 0{N'^) just as for the 
local dimer flip dynamics treated in the present paper. 

The height fleld is the hidden order parameter in the 
present model (or of height models in general) and hence 
is the proper way to approach the coarse-grained, long 
time dynamics. Note that experience has shown that 
the static exponents are extracted much more efficiently 
from (|/i(q)p) than from spin-spin correlations (using the 
same simulations) 1^,^,^. Therefore, I suggest numerical 
studies of the dynamics ought to extract estimates of the 
correlation times directly from the height fields rather 
than from spin-spin correlations. 
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APPENDIX A: CORRELATIONS OF SPIN 
OPERATORS 

This appendix shows that a standard order parameter 
m{r) of dimer coverings [p|j2l|] reduces to a function of 
the heights, and correspondingly its asymptotic correla- 
tions are a function of those of the heights. (This works 
for the standard order parameters of all height models; 
the derivation for equal-time correlations is discussed at 
greater length in Refs. ll]j|-|fl.) 

Each site r = (x,?/) of the direct (not dual lattice), is 
connected to one other site r' — {x',y') by a dimer., De- 
fine m(r) ~ (—1)^+^ +1 ((_i)J'+y +1^) for dimers running 



^It turns out that all the algorithms cited have a simple 
action on the height variables, thus the existence of a height 
map may be a prerequisite to the possibility of accelerating 
such models. 
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in the ±x (±?/) direction.|^ Then (with the arbitrary 
constant offset of z(r) fixed as in Fig. |l|) 



(Al) 



where h(r, t), defined on direct sites, is the average of the 
four z{r,t) values on the surrounding dual sites. 
Then 



G(r,t)= (TO(0,0)m(r,i))«exp[-i(^-j C{r,t)] (A2) 
where 



(A3) 



C{r,t)=([hir,t)-hiO,t)f) 



Using Fourier transformation with eqs. (|lj) and (|13|), one 

obtains 

= / £ " e"'^'"''''* cos(q • r)], (A4) 

a logarithmic divergence cut off at high q by the inverse 
lattice constant and at small q by min(r~^, [TKt]~^^'^. 
Thus finally the correlation function behaves as 



1 



G{r,t) w —^{r/VVKi) 



(A5) 



' tt/2 



where $() is a scaling function, such that G{r, t) ^ r ''"/^ 
or as G(r, t) ^ (riirt)"^''-/^ depending whether r/VVKt 
is much less than or much greater than unity. This shows 
that the dynamic critical exponent defined from the stan- 
dard correlation functions is indeed z = 2. |^ 

APPENDIX B: STATIC FLUCTUATIONS 

This appendix reviews what is know exactly about the 
equal-time fluctuations of Fourier components of z(r). 
This result will be an essential lemma for the proof about 
A(q) outhned in Sec. 0. 

Evidently, in the continuum picture, the equal-time 
expectations are given by 



q,q' 



K\q\'' 



(Bl) 



(Recall that h* = h^q.) This scaling is equivalent to 



{\h{r) - h{r')\') = 0(1) + ^ In |r - r'| (B2) 



One could equivalently consider m(r) as living on the 
occupied bonds, since m(r) always takes the same value at 
either end of a dimer. 
^^Eq. A5 was previously derived in the context of the hexatic 



order parameter: see eq. (5.13) of Ref. 



in direct space. 

In the study of "height models" [^,||j|J^ and random- 
tiling quasicrystals p7[| , it is the fundamental hypothesis 



that i(q) has a behavior like (Bl) at small q. The be- 
havior (Bl) has been proven rigorously for a very lim- 
ited set of discrete models Q, and not even for the 
well-known BCSOS model of rough interfaces pO| . Even 
though the latter model is "exactly soluble" [[4l[ via the 
Bethe ansatz, such correlation functions are not known 
exactly (correlation functions are notoriously difficult to 
extract from the Bethe ansatz method p2)). 

However, the dimer covering model (with 2D Ising 
models) belongs to the "free-fermion" class of exactly 
soluble models |H,|l^, and using the Pfaffian approach 
the dimer-dimer correlation functions can (in principle) 
be written out exactly Later on McCoy evaluated 

the large-distance asymptotic behavior of the correlation 
function, for any orientation of the vector between the 
two dimcrs 

Say that Ux (r) = 1 when there is a dimer connecting r 
tor+[l,0] and zero otherwise; ^^^^(r) = nx{r)—{nx{r)) = 
'T-x(r) — 3, similarly for ny{r). 

It was computed [|l5| that, in the equal- weighted en- 
semble of dimer coverings, 



CxxiR) = {Snx{r)Snx{r')) = 
-g{X,Y)^+g{X + l,Y)g{X -1,Y) 



(B3) 



Here R = {X, Y) = r' — r, and gjX, Y) is a Green's 
function arising from the Pfaffians At large R, 



- — X odd, Y even; 
g{X,Y)-{iX^^ X even, r odd; (B4) 




X even, Y odd; 
otherwise. 



Hence, collecting the four cases {X and Y even or odd), 

1 , 



-1 



(B5) 



Dimer correlations imply correlations in height gradi- 
ents since the definition of z(r) (Sec. 0) is equivalent to 



^(r+[|,i])-^(r + 
^(i-+[|,i])-^(r + 



i~lf+HSny{r) (B6a) 
{-iy+y+HSnx{r) (B6b) 



Combining (|B5[) and (B6b), with 



z(r+[i,i])-z(r+[i,-i])^ V,z(r) 



(B7) 



gives 
(V,z(r)V,,z(r')) 



The term in (B8a) multiplied by (—1) will be neglected 
since it averages to zero. Similarly, using the analogous 
formulas for {SuxSuy), 
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42 2XY 
27r2 i?4 



(B8b) 



and 



(V.z(r)V..z(r')) = 



(B8c) 



The non-oscillating correlations in Eqs. ( |B8D have the 
form of dipole-dipole interactions, as if Vz(r) is a polar- 
ization. | p3[ Integrating these equations with respect to 
R, 



16 



(|z(r) - z(r')P) ~ const +^ In |r-r'| 

and Fourier transforming yields the required result, 

16 



(B9) 



7r|q|- 



(BIO) 



The behavior is indeed described by the continuum mode 
@, with K = 7r/16. 

An exact result for arrow-arrow correlations was also 
obtained for a special case of the 6- vertex model p4[ , 
and applied to compute exactly the coefficient of the log- 
arithmic asymptotic behavior of (|z(r) — z(r')p) for the 
corresponding BCSOS model @. (The BCSOS model is 
just the height mapping of the 6- vertex model). Actually, 
the parameter values in that special case make it equiva- 
lent to a random dimer covering. To convert the h{r) of 
the dimer model to the height hBcsosi^) of ^ BCSOS 
model (with lattice constant V^), put h{r)/2 — hscsos 
on the sublattice with even ft,(r) and erase h{r) on the 
sublattice where it is odd. (See also [|5|; to relate this to 
other versions of the same mapping between models ||4^ ] 
reverse all the arrows pointing along vertical bonds). 
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FIG. 1. (a). Random dimer packing, showing heights z{r) 
(those in flippable plaquettes are circled), (b). Elementary 
dimer flip, showing heights. 
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